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Abstract Nonequilibrium dynamics of an A/'-fold 
spin-degenerate ultracold Fermi gas is described in terms 
of beyond-mean-field Kadanoff-Baym equations for cor- 
relation functions. Using a nonperturbative expansion in 
powers of 1 /Af, the equations are derived from the two- 
particle irreducible effective action in Schwinger-Keldysh 
formulation. The definition of the nonperturbative ap- 
proximation on the level of the effective action ensures 
vital conservation laws as, e. g., for the total energy and 
particle number. As an example, the long-time evolu- 
tion of a homogeneous, twofold spin-degenerate Fermi 
gas is studied in one spatial dimension after an initial 
preparation far from thermal equilibrium. Analysis of 
the fluctuation-dissipation relation shows that, at low 
energies, the gas does not thermalise. 
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C\| 1 Introduction 



■ The preparation of ultracold atomic Bose and Fermi 
' gases in various trapping environments allows to study 

■ important aspects of quantum many-body dynamics. 
Many of the past experiments with ultracold gases can 
be approximately described by semi-classical approxima- 
tions of the quantum many-body field equations such 
as by the Gross-Pitaevskii, Hartree-Fock-Bogoliubov, or 
Bardeen-Cooper-Schrieffer theories [Tl[2]. The descrip- 
tion of the dynamics of many-body systems of suffi- 
ciently weakly interacting particles is usually based on 
perturbative approximations that rely on an expansion 
in powers of some dimensionless parameter that mea- 
sures the binary interaction strength. These approxima- 
tions are generically based on the smallness of statistical 
fluctuations. In the limit of infinitely strong interactions, 
a number of approaches exists (in particular for systems 
in only one or two spatial dimensions) that allow dual 



descriptions in which approximations rely on the small- 
ness of the inverse of the coupling strength [3111]. 

For intermediate interaction strengths, only a few ap- 
proaches exist. For such systems, quantum as well as 
strong classical fiuctuations play in general an impor- 
tant role. Prime examples are given by ultracold gases 
driven into the vicinity of Feshbach-Fano scattering res- 
onances [5j|6j[7l[8] , lattice- trapped gases in between the 
superfiuid and Mott insulator regimes [Ql fTOlfTTj , and low- 
dimensional gases in regimes where the transverse con- 
finement strongly affects the binary scattering dynam- 
ics of the particles pi [T^[T5] . Ultracold Fermi gases have 
been studied extensively in the vicinity of the BEC-BCS, 
i.e., superfluid-superconducting crossover [l4l[T5lfT6lfT7l 
[T51ll9| . and currently attract increasing interest, e.g., 
in the context of Kondo phenomena in lattice environ- 
ments [2D] . 

Taking into account higher-order classical and quan- 
tum fluctuations is important in describing the late- 
time behaviour of initially strongly perturbed as well 
as of continuously driven systems. A quantitative de- 
scription of the field dynamics in a closed system 
then generically requires nonperturbative approxima- 
tions that need to be consistent with vital conservation 
laws such as the conservation of the total energy and 
the total particle number. In this article, we present a 
self-consistent formulation of the nonequilibrium dynam- 
ics of an A/'-component ultracold Fermi gas in terms of 
beyond-mean-field Kadanoff-Baym equations for corre- 
lation functions. Dynamic equations for the two-point 
correlation functions are derived from the two-particle 
irreducible (2PI) effective action [2n[22ll23] in Schwinger- 
Keldysh formulation [241125] . This approach extends the 
self-consistent mean-field formulation by including the 
effects of scattering between particles to the desired or- 
der of approximation and ensures the local conservation 
of the particle density-current vector as well as of the 
energy-momentum tensor. Self-consistency is achieved 
by invoking approximations on the level of an effective 
action and obtaining time evolution equations by subse- 
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quent functional derivation with respect to the desired 
correlation functions. 

Functional integral techniques based on the 2PI ef- 
fective action enable an efficient embedding of sum- 
mations of infinite series of perturbative processes. In 
recent years, systematic nonperturbative expansions of 
the 2PI effective action to next-to-leading order in in- 
verse powers of the number of field components I26| 
have allowed substantial progress [27,.28. . After first 
successful applications of these nonperturbative expan- 
sions to the study of far-from-equilibrium dynamics as 
well as thermalization in relativistic bosonic ^27.,29,.30, 
and fermionic [311132) theories, they have recently been 
employed in the context of ultracold bosonic quan- 
tum gases [33l[34l[35ll36ll37ll38] . A related alternative ap- 
proach based on renormalization-group-like fiow equa- 
tions was presented in Ref. [5Ull40j . For introductory 
texts see, e.g., Refs. [ITlli^ . 

Considering the case of a degeneracy in the M spin 
degrees of freedom, we derive the 2PI effective action 
to next-to-leading order (NLO) in a nonperturbative ex- 
pansion in powers of l/M . While the l/M approximation 
is entirely classical in leading order, classical fluctuations 
and corrections induced by quantum fluctuations, both 
of which are introduced by scattering, are included in 
NLO in a nonperturbative manner. A similar l/M ex- 
pansion can be derived for SU (Af) symmetric models of 
Kondo lattice systems |20) . 

After introducing our approach in detail, we present 
in this paper, as an example, the long-time evolution of a 
homogeneous, twofold spin degenerate Fermi gas in one 
spatial dimension following a quench in the interaction 
strength. Our results extend upon the work presented in 
[43] . The considered model is integrable and has as many 
conserved quantities as there are degrees of freedom [44] . 
Hence, it is expected not to thermalise in general [151H51 
113. 

Its low-energy properties can be approximated by 
a Tomonaga-Luttinger liquid (TLL) model 05] which 
contains a linear free dispersion. The resulting model 
is known to form a low-energy flxed point of the full 
interacting one-dimensional Fermi gas. Owing to the 
quadratic form into which the TLL Hamiltonian can be 
transformed by introducing bosonic particle-hole oper- 
ators, the occupation numbers of the resulting quasi- 
particle modes represent conserved quantities. In A9, 
154) . the long-time evolution of TLL fermion and cou- 
pled fermion-boson models after an interaction quench 
was found to approach a generalised Gibbs ensemble [45] , 
accounting for the conserved quasiparticle numbers. We 
point out that non-thermal stationary states have been 
found also in a number of other models, both integrable 
and non-integrable [50l[5n[52ll53) . 

Here, we consider the dynamic evolution described 
by the full interacting fermionic Hamiltonian, approach- 
ing the problem of equilibration from the high-energy 
end. The KadanofF-Baym dynamic equations in 2PI NLO 



l/A/" approximation we employ are generically consid- 
ered to describe thermalisation I27ll29[ l30ll3n[32l[33ll34l 
[35l[36ll37] . While the conserved quantities of the consid- 
ered ID Fermi gas are not known explicitly, and while we 
do not expect to recover the TLL model at very low en- 
ergies, we flnd, by analysing the fluctuation-dissipation 
theorem, that the correlation functions emerging at late 
times are incompatible with a thermal ensemble. 

Our paper is organised as follows: In Sect. [21 we 
deflne the fermionic model Lagrangian and summarise 
the most relevant aspects of the 2PI effective action ap- 
proach to nonequilibrium dynamics. Integro-differential 
Kadanoff-Baym or Schwinger-Dyson dynamic equations 
for the two-point Green functions are obtained from the 
effective action. In Sect. [3J we present different non- 
perturbative approximations of the 2PI effective action 
based both on a coupling expansion and on an expan- 
sion in powers of l/A/". As an exemplary application 
demonstrating the power of the method, we study, in 
Sect. 21 the long-time evolution of a homogeneous one- 
dimensional Fermi gas with twofold spin degeneracy 
from an inital state that is far from thermal equilibrium. 
Our conclusions are drawn in Sect. |S] In the appendices, 
we review relevant properties of Grafimann variables and 
two-time Green functions, and provide a discussion of 
the conservation of total energy and particle number by 
the approximated 2PI effective action. 

2 Two-particle irreducible effective action 
approach to nonequilibrium dynamics 

2.1 The model Lagrangian 

We consider the dynamical evolution of an ultracold 
Fermi gas of atoms far from thermal equilibrium. The 
atoms are assumed to be, internally, in Af different hy- 
perfine states. Hence, in the language of quantum held 
theory, we study the dynamics of nonrelativistic complex 
fermionic flelds !Z'q(x, t) obeying canonical anticommuta- 
tion relations 

[#„(x,t),i^^(y,i)]+ =<5„/3<5(x-y), (1) 

where [•, •]+ denotes the anticommutator. The indices 
a and /? count the Af hyperfine "spin" states, s-wave 
contact interactions between atoms in different hyper- 
fine states are assumed, while Pauli's principle forbids 
s-wave collisions between fermions that are internally 
in the same state, p-wave and higher-order partial-wave 
contributions are neglected. Our formalism, though, can 
be readily extended to more complicated interaction po- 
tentials. In the contact potential approximation, the in- 
teractions in the channel characterised by the asymp- 
totic hyperfine states a and /3 are described by the po- 
tential T^nt,a/3(x - y,t) = ga^(i)(5(x - y) with a pos- 
sibly time-dependent coupling strength gafs{t). In three 
spatial dimensions, the coupling strenght is related to 
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the scattering length between states a and /? by 
the relationship gap = 'iT^a-ap /m. Hence, the interaction 
Hamiltonian reads 

H,^,{t) = i^ljli^)rpl(^)4rp(^)^a{^), (2) 

where J^ = J d'^x in d spatial dimensions. 

In this article, we will be concerned with a func- 
tional integral formulation of the dynamical field the- 
ory that involves functional integrations over complex 
Grafimann- valued fields ipaix) obeying 



imaginary parts of the quantum field. 



(5) 
(6) 



To simplify the notation, we include the hyperfine index 
a and the field index i into a single index a = (a, i). 
Sums over a imply a sum over a £ {!,•■• , A/"} and one 
over i e {1,2}. 

The Grafimann action 



(7) 



[ijaix),^piy)]+ = [i^a{x),ip}{y)]+ = (3) ^^^^ = J dxo, associated with the Lagrangian Q 



for any combination of a, /3, x, and y, instead of a formu- 
lation in terms of equal-time anticommuting field oper- 
ators. Here, x = (t,x) = (a;o,x), etc., denote space-time 
coordinates and the asterisk complex conjugation. Prop- 
erties of Grafimann variables relevant for our discussion 
can be found in App. 

The Lagrangian for the A/'-component ultracold 
Fermi gas with the above interactions reads 



1p*aix)dxo'lpa{x) - [d^„1pa{x)]tl;a{x) 



(4) 



+ ^*a{x)^i^aix) - 1pa{x)Vc^t,a0{x)^pp{x) 
- ^^aixWpix)i:p{x)^a{x] 

where dxg denotes the partial derivative with respect to 
time, and Vcxt,ai3ix) are possibly time-dependent trap- 
ping potentials or other external-field matrix elements 
coupling the hyperfine levels a and /3. Summations over 
repeated indices are implied. A factor of 1/Af has been 
taken out of the couplings Xap = Mgap in order to make 
the relative weight of the interaction to the quadratic 
terms in the Lagrangian invariant under a rescaling of 
M . This will be of use when considering the expansion 
in powers of \/N in Sect. 13.21 

Three symmetries play an important role: Local con- 
servation of particle number implies that the Lagrangian 
density possesses a global C/(l) symmetry within each 
hyperfine subspace. Moreoever, if the couplings between 
the hyperfine levels through both the external field 
^4xt,a^ and the interactions Xap are all equal to each 
other symmetry under rotations in the space of hyper- 
fine states. Furthermore, for vanishing external fields 
and time-independent interactions, the Lagrangian is 
Galilei invariant, implying a locally conserved energy- 
momentum tensor. 

For the ease of numerical implementation, we choose 
a field basis "^q.i where the index i distinguishes real and 



reads 



1 



= 7^ / ^a{x)iG-^,,{x, y)My) + ^intM • (8) 

^ Jxy 

Here, the inverse free fermionic propagator is given by 

'^Go,lbix,y) = S{x~ y){\Tabdxa - KttH^Jpix)) , (9) 
with the one-body Hamiltonian 
V2 



H^Jlix) 



2m 



Sap + Vcxt,ap{x) 



(10) 



Here, J = J dxo J d'^x denotes the integration over the 
region of space-time under consideration. Sab — SapSi^i^, 
and 



_ r 2 2 _ /O -i 

Tab = -Sap(Ti^i^ , ^ ^ I i 



(11) 



Furthermore, S{x — y) = 5{xq — yQ)S^'^\yi — y) denotes 
the (d + l)-dimensional Dirac distribution and ipaix) = 
4'b{x)Tba- The interaction part Sinti'ip] corresponding to 
the Lagrangian Q is 



Sint bP] = - 



XaP 

8Af 



1pa{x)lpa{x)-)pbix)'>Pbix) . (12) 



Note that for a Fermi gas trapped in a lattice po- 
tential, the action, in the tight binding approximation, 
has the same form as in Eq. ([5]), and the free inverse 
propagator is given by 



iG'o,a6(2;> y) = iTabd.x,,S{x - y) 

- SabHa^{x, y)S{xo - yo) , 



where 



H^''{x,y)^~JS, 



id) 

(n,m) 



(13) 



(14) 



^ If not otherwise stated, we use natural units, with h= 1. 



is the one-body Hamiltonian. Here, x — (xo,n) and 
y — (yojin) denote the lattice space-time coordinates, 
S{x-y) = SnmS{xQ-yo), and s\^'^^^^ = 1 if and only if, in 
the single-band approximation, n and m denote adjacent 
sites in the d-dimensional lattice; otherwise, <^|n^rn) ~ 
The site dependent energy e^.n describes, e. g., an addi- 
tional external trapping potential for hyperfine mode a, 



4 



Matthias Kronenwett and Thomas Gasenzer 



and we neglect spin mixing by this potential. All other 
previous and subsequent equations carry over to the lat- 
tice case when spatial integrals are replaced by sums over 
the lattice sites. Extensions to the case of two or more 
internal electronic states and more specific symmetries 
in the hyperfine levels, as used, e.g., to describe Kondo 
lattice systems [5D], are straightforward and will be the 
subject of a subsequent publication. 

2.2 Nonequilibrium generating Junctional 

Knowing the time dependence of the correlation func- 
tions of a many-body system allows to derive the dy- 
namics of physical observables. These correlation func- 
tions can be obtained from the nonequilibrium generat- 
ing functional Z[K\pi:)\, 



PD{to)Tce' 



(15) 

where poita) is the normalised density matrix of the 

system at the initial time to, and ^a{x) — tf'f,(x)Tf,a, 
^a{x) are field operators in the Heisenberg picture. 
C = C'^ U C~ indicates that the temporal integrals are 
taken along the closed (Schwinger-Keldysh) time path 
(CTP) [241125] from the initial time to to infinity (path 
C"*") and back to to such that ^ — Jq and 

Jq ~ J(2+ — dxQ. Tc denotes time-ordering 
along the CTP, which implies that operators evalu- 
ated at later times stand to the left of those evaluated 
at earlier times. The classical external two-point field 
Kab{x,y) = ncKcd{y,x)Tda is introduced to allow for 
the generation of correlation functions of order 2n, 



{Tcl'a, (yi) • • • <f a„ iXn)^b^ iVn)) 

2" 5^Z[K;pD{to)] 



Z i5Ka^bt{xi,yi) ■ ■ ■i5Ka^bSxn,yn) 



(16) 



K=0 



where (•) = Tr[pD-]. Translated into a functional inte- 
gral, the generating functional takes the form 

z[K-pD{to)] = / v^+v% {4>+\pD{tQ)\4>^) 



i^a{x)Kab{x,y)lpb{y) 



C,xy 



(17) 



Here, \4>^) are eigenstates of binary field operator 
products evaluated at the beginning and the end 
of the CTP, respectively, ^„(x, to)#a(x, to)|0^j=) = 
-i^^(x,to)TbaiA^(x,to)|0j)- The differential 
parts of the path integral measures are defined 
as = na,xdV'a,i(^'^o)dV'a,2(x,io) and 

= na,x,a;oec\{to}'^^".i(^)'^^«.2(a;). The action 
on the CTP is defined as 



C,xo 



(18) 



The most general initial density matrix can be 
parametrised as [4ni42| 

i^^^lMtoMlr^) =AAoexp[i/c[^]] (19) 

with normalization factor Aq and a functional fc ["0] that 
can be expanded in powers of the fields, 

/cN=a("^+EA / ait.aAx\...,xn 

n 

X n^a,^^:-), (20) 

m— 1 

where the above boundary conditions are implied. The 
cumulants ai"La^(a;^, . . . , a;") are non-zero only at time 
to where the density matrix poito) is specified. In the 
following, we will only consider Gaussian initial states, 
for which the cumulants vanish for all n > 2 also at time 
to. This allows to combine the initial density matrix with 
the external source field by defining the nonlocal source 

^abi^^y) = Kab(.^'y)-Tacacb {^,y)- This allows to Write 
the generating functional in the simpler form 



Z[K] = J Vi^ expji^cM 



i'a{x)Rab{x,y)i^b{y) 



Cxy 



(21) 



where the measure also includes the fields at time to. 

The connected two-point Green functior|^ in the pres- 
ence of the nonlocal source R can be derived by func- 
tional differentiation, 



-Gba{y,x\R) = —— 

2 dRab{x,yj 



of the Schwinger functional 



W[R] = -ilnZ[i?] 



(22) 



(23) 



which is the generating functional for G. Later on, we 
use a vanishing external source K = 0, and employ the 
notation, cf. Eq. p^ . 

Gab{x,y) = (TcMxWbiy)) = Gabix,y;R[K = 0]) . 

(24) 

Symmetry properties of the two-point function G as well 
as its decomposition into the physically relevant statis- 
tical and spectral correlation functions are discussed in 
detail in App. |B] 



^ For fermionic fields, the two-point Green function G is 
identical to the connected two-point function since the field 
expectation value always vanishes. For higher n-point func- 
tions, there is, however, a distinction. 
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2. 3 2PI effective action 

Directly evaluating the full quantum real-time path in- 
tegral in Eq. ((2T|) is in general not feasible as the oscil- 
lating complex measure represents a variant of the sign 
problem. One therefore needs to resort to analytical ap- 
proaches in evaluating the dynamics in regimes where 
quantum fluctuations are relevanl|f|. 

The goal in deriving an effective action is similar to 
that of classical mechanics in the Lagrangian formula- 
tion: The effective action r[G\ is defined such that it 
allows, by Hamilton's principle, to derive the dynamic 
equation for the correlation function G. The action takes 
into account quantum effects, and the dynamic equations 
derived from it will obey crucial conservation laws cor- 
responding to the symmetries of the underlying model 
Lagrangian. 

The two-particle irreducible (2PI) effective ac- 
tion pni22l[^ for the fermionic Lagrangian ^ is de- 
fined by a Legendre transform of W[R\ with respect to 
the source 



r[G] = W[R] - [ 

Jxy 



SW[R] 
5Rab{x,y) 



Rab{x,y) 



^w[r] + -i:y[gr] , 



(25) 



where we have used Eq. (j22p . and where it is implied 
that Eq. (l22t can be inverted to give i? as a function of 
G. The effective action r\G] can be written as a series 
of terms represented by two-particle irreducible vacuum 
diagrams [3T], 

r[G] =-^Tr[ln(G-i)+Go^G]+/^2[G]+const.. (26) 

The first term is a one-loop-type term derived from a 
saddle-point approximation of the path integral, where 
the trace denotes summation over all field and spin in- 
dices, and integration over spatial coordinates and over 
times along the CTP. The constant term is irrelevant for 
the dynamics. While the one-loop term remains within 
the mean-field approximation, scattering effects are con- 
tained in r2[G]. 12 [G] can be written as the series of 
all vacuum, i.e., closed, 2PI diagrams constructed from 
the Green function G and the bare vertices defined by 
the model Lagrangian pi]. A diagram is 2PI if it does 
not fall apart on opening two of its lines. The specific 
expansion used in this work will be described in Sect. [3l 



^ This is in particular the case for long-time evolutions and 
where interactions become strong. If quantum fluctuations 
are small then the quantum part of the fluctuating fields can 
be integrated out leading to a classical path integral that can 
be computed using Monte Carlo techniques [3611551156 II57II58) . 

Note that, as compared to the Bose case, there is no Leg- 
endre transform with respect to the one-point source Ja{x) 
as there is no field expectation value, either. 



2.4 Dynamic equation for the Green function 

Given the effective action F, the stationarity condition 

6r[G] 



5GUy,x;R) 2 



-Rab{x,y) 



(27) 



determines G for a given nonlocal source R. For a given 
initial state, Eq. (P7|) represents the equation of mo- 
tion for the Green function G. Using Eq. (pS)) . one finds 
the real-time Schwinger-Dyson- or Kadanoff-Baym-type 
equation 

Gab{x,y\R) = G~l^{x,y) - iRab{x,y) - Eab{x,y;G) , 

(28) 



where 



i:abix,y;G) = -2i- 



<^r2[G] 



dGha{y,x;R) 

denotes the one-particle irreducible (IPI) self-energy. 
Convolving Eq. (1^51) with G yields the dynamic equa- 
tion for the two-point Green function, 

Gq\c{^, z)Gcb{z, y; R) = Sc{x - y)Sab 
+ [ {S,c{x,z;G)+iRacix,z))Gcbiz,y;R). (30) 

Jc,z 

For a closed system (i. e., vanishing external source K = 
0), no explicit dependence on R remains in the evolution 
equations since the term ^ iRadx, z)Gcb{z, y; R) van- 
ishes for xq ^ to, i.e., the term becomes redundant for 
all times of interest. 

Note that the free inverse propagator, Eq. con- 
tains a first-order time derivative. Since the free inverse 
propagator is otherwise diagonal in (x — y), the inte- 
gral on the left-hand side can be carried out, revealing 
the integro-differential structure of the dynamic equa- 
tion ((30t . Despite Eq. ((30|) being exact, to be solved, 
it requires knowledge of the self-energy and therefore of 
the 2PI part /^^ of the effective action. For practical com- 
putations, truncations of the series of 2PI diagrams are 
chosen as discussed in the following sections. 



3 Nonperturbative approximations of the 2PI 
effective action 

In this section, we discuss different possible truncations 
of the expansion of / 2 in terms of two-particle irreducible 
(2PI) diagrams. Figure [T] shows the leading diagrams in 
the series ordered by the number of bare vertices per 
diagram. The Green functions are represented by (blue) 
solid lines, the vertices by (black) dots. 

For the following discussion, it is convenient to sepa- 
rate out the local contributions to the proper self-energy. 



Eab{x,y;G) = -iZ'^"\a;; G)(5c(a; - y) 
+ Sab{x,y;G) , 



(31) 
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r2[G\ 




+ 



+ 



^ 2 



[G] 




Fig. 1 (Color online) Diagrammatic expansion of the two- 
particle irreducible (2PI) part i~2 of the effective action (|26p 
in terms of 2PI graphs. (Blue) solid lines stand for the Green 
function G, and (black) dots for the interaction vertex A. 
Explicitly shown are all diagrams that contain up to four 
vertices. All statistical factors and relative minus signs are 
omitted. 



a 



a 





+ 




+ 




a 



(3 p [3 p 



a 



Fig. 2 (Color online) Decomposition of the bare vertex A 
(black dot) as a sum of the three possible spin-index contrac- 
tions. The point interaction vertex Xap is represented by a 
squiggly line at each end of which spin and field indices are 
conserved and summed over. 

and include them together with the one-body Hamilto- 
nian term of G^^ in the matrix 



Mab{x,y-G)=5c{x-y)[5,^,,H'J^,{x) 



(32) 



The dynamic equation (pO| . for if = 0, now takes the 
compact form 

iTacdxaGcb{x,y) - i6abSc{x ~ y) 

{Mac{x, z- G) + iSac(x, z; G)) Gefc(z, y) , (33) 

where Gab{x,y) is defined in Eq. Note that, ne- 

glecting S, the dynamic equation reduces to a local dif- 
ferential equation in time, while the additional nonlo- 
cal contribution to the integrand introduces memory of 
the past evolution of G into the equation. Furthermore, 
it is convenient to reveal the specific structure of the 
point interaction vertex and to represent the vertex by 
a squiggly line. The interaction part of the considered 
Lagrangian density ^ requires that spins are conserved 
at each end of the vertex. Furthermore, it does not allow 
the Green functions in two adjacent loops to correspond 
to the same hyperfine quantum number since Xaa = 0. 
The three possible index contractions at each vertex are 
depicted in Fig. [21 



3.1 Coupling expansion 

The expansion of I2 in terms of 2PI diagrams can be 
ordered as a power series in the bare coupling constant 



^2nd-order j'^j 



pHFB ^ 




+ 2 




Fig. 3 (Color online) The first-order (Hartree-Fock- 
Bogoliubov) and second-order coupling expansions of /2 • The 
factors of 2 arise due to the different possible spin index con- 
tractions depicted in Fig. (2] All other statistical factors and 
relative minus signs appearing because of the Grafimann al- 
gebra properties are omitted. 



g. Note, that this, in the strict sense, does not consti- 
tute a perturbative expansion in powers of g since the 
Green function G entering the diagrams is determined 
self-consistently by Eq. ([50)1 . Thus, each diagram con- 
tains contributions up to arbitrarily high powers in the 
coupling. Nevertheless, if the solution of the dynamic 
equations exhibits G to be a bounded function then the 
coupling expansion of /2 can be effectively viewed to be 
perturbative. 

3.1.1 Hartree-Fock-Bogoliubov approximation As a 
first step, we recover the dynamic equations in the 
Hartree-Fock-Bogoliubov (HFB) approximation. Re- 
taining only the lowest-order diagram of the coupling 
expansion of /2, that is, the double-bubble contribution 
shown in Fig. [3l is known as the Hartree-Fock- 
Bogoliubov approximation. Due to the structure of 
the vertex, Eq. p^ . this contribution consists of two 
qualitatively different diagrams and reads 



Aq/3 



Gaa{x,x)Gbbix,x) 



2Gabix,x)Gba{x,x) 



(34) 



where it is summed over a = (a, ia) and h = (/3, if,). The 
self-energy derived from Eq. ([M)) is 

E^rix, y; G) ^ ~iS^r^'\x; G)5c{x - y) , (35) 

ET''\--^G)^~^(5.bGU-,x) 

2-^ ^ (36) 
- 25^iiGab{x,x) \ , 



where it is summed over c = (7, ic)- Hence, the HFB 
self-energy reduces to its local part £'^^^(0). 
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We employ the decomposition of the fuh Green func- 
tion G into the statistical correlation function F and the 
spectral function p, 



Gabix, y) = Fab{x, y) - -pab{x, y) sgnc(a;o - yo) , (37) 



which are defined in Eqs. (|108p and (|109l) . respectively, 
and discussed in more detail in App. |B1 Inserting this 
decomposition into the dynamic equation for G and us- 
ing Eq. (I113P , we find the time-dependent Hartree-Fock- 
Bogoliubov equations for p and 



\''^Tac^xo 



M, 



HFB 



[^Tacdxo - Mac 



HFB 



(a:;G)]F,b(x,y) = 0, (38) 
(a;;G)]pefc(x,2/) = 0, (39) 



where Af"^^(a:;G) is obtained by Eq. with 
G) = S^^^^°\x; G) as given in Eq. As the 
nonlocal self-energy contribution vanishes in the HFB 
approximation, the dynamic equations for F and p de- 
couple. At equal times, xq = j/Qj the spectral func- 
tion p is fixed by the anticommutation relations (jll2p . 
Hence, one finds that the single-particle density matrix 
nai3{x,y,t) ~ (!f'^(x)!f'^(y))f is solely determined by the 
statistical correlation function: 



"q/3(x, y, t) - -SapSix - y) 

- -i<5,„i5,,i(Faf,(x,y,t)+^^ab(x,y,i) (40) 
+ i{F^-,iK,y,t) - Fab{x,y,t))^ , 



where a — (a, 3 — ia), and it is summed over ia and ib- 
This includes the density of particles in mode a at point 
(x, i), na(x,t) = nQQ,(x, X, t). Moreover, the anoma- 
lous density matrix or pair function r7iQ,^(x, y, t) = 
(#Q(x)#/3(y))t is given as 



TOa/3(x,y,t) 

= -^3iai^i,i[Fab{x,y,t) - Fg^-b{x,y,t) (41) 
+ i{F^bix,y,t)+Fab{K,y,t))) , 



where it is summed over ia and if,. Adding Eq. (1381) and 
its transpose, one obtains the set of coupled HFB equa- 



tions for nQ^(x,y, i) and mQ^(x,y,i): 

{idt + Haaix) - Hpi3{y))nai3{yi,y,t) 

= 1 ^;^(~%7(x,x,i)na^(x,y,i) 

-|-TO*^(x,x,t)m^/3(x,y,t) (42) 

-f- fia^f (x, X, t)n^f3 (x, y , t) 

- {(a,x) O (^,y)}* , 
{idt ~ Haa{x) - Hpi3{y))map{-x,y,t) 

(n^j (x, X, t)maf3 (x, y, t) 

- n*a-y (x, X, t)m^p (x, y , t) (43) 

- TO„^ (x, X, t)h-^p (x, y, t) 



-{(a,x)^(/3,y)}, 

where it is summed over 7, and nap{'x.,y,t) = 
na^(x, y,t) — (5q,^J(x — y)/2. The last term in Eq. (03]) 
(Eq. (HH)) denotes (the complex conjugate of) the first 
term in curly brackets with a and /3, and x and y in- 
terchanged. Equivalently, the HFB equations can be de- 
rived using the Ehrenfest theorem, i. e., i9tf^a/3(x, y, t) = 
-([iJ, tf't(x)tf'^(y)]_) and similarly for ma/3(x,y,t), and 
a Weyl ordered Hamiltonian H. When evaluating the ex- 
pectation values, the Hartree-Fock-Bogoliubov approxi- 
mation then consists of neglecting all joint cumulants 
higher than second ordei[l. 

In this paper we will consider, as a concrete example, 
a homogeneous one-dimensional Fermi gas "on a ring" , 
i.e., in a finite box with periodic boundary conditions. 
In this case, the correlation functions only depend on the 
spatial relative coordinate x — y. In this case, Eqs. ([55)) 
and pop are conveniently solved in momentum space. 
The solutions are 

Fab{v,t) = (exphirAfH™(p)t])^^i^,,(p,0 , (44) 
Pab{v,t) = (exphiTAf"FB(p)i])^^^^^(p^i) ^ (45) 



with 



Mfr{P)=5ab 



2m 



(46) 



where = (27r) J and it is summed over 

c = (7,ic). Inserting the solution for F into na{p,t) = 
[1 - -f^(Q,i)(a,i)(P' - P{a,2){a,2){p, t)]/2, wc rccovcr that 
the HFB equations, which exclude scattering, leave all 
momentum-mode occupation numbers invariant. 



^ For related discussions in the context of cold gases, see, 
e.g., Refs. [34||59h60,|61(|62]. 
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3.1.2 Second-order coupling approximation Beyond the 
mean- field Hartree-Fock-Bogoliubov contribution to /2 , 
we now additionally take into account the second-order 
diagram shown in Fig. [3] containing two bare couplings: 



with 



Gbc{x,y)Gcbiy,x) 



xy 



X Gda{x,y)Gadiy,x) (47) 



j Gad{x,y)Gdb{y,x) 



X Gbc{x,y)Gca{y,x) , 

where it is summed over a = (a,ia), b = {/3,ib), c = 
(7,ic), and d — {S,id). Taking the functional derivative 
with respect to G, the additional term yields the nonlocal 
self-energy 



^abix,y) = '^oIr2^ n^six, y)Gab{x, y) 



2AA2 



(48) 



Gad{x, y)Gdciy, x)Gcbix, y) 



with 



nal3{x,y) = G(^a,i)if3,j){x,y)G(^f3j)(^a,i){y,x) , (49) 



where sums over 7, S, ic, and id are implied in Eq. 
and over i and j in Eq. (HH)) . 

In order to derive the contributions to the dynamic 
equations for F and p, one decomposes S and 77 into 
their statistical and spectral parts, 



Sab{x,y) = Sab{^,y) - ^^^f,(a;,y)sgnc(a;o -yo) , 



(50) 



n„l3{x,y) = n^fi{x,y) - -n^i^{x,y)sgnc{xo ~ yo) 



which gives 



^abi^^y) 



2AA2 



(51) 



n^si^,y)Fab{x,y) 

-\nPg{x,y)pab{x,y) ^^2^ 
-2P[scix,y)Fcbix,y) 

+ 2^aSci^^y)Pcbix,y) 



V.P .A ' ^°]^2^ ( n^s{x,y)Pab{x,y) 



+ nP,{x,y)Fab{x,y) 
- 2P[sc{x,y)Pcb{x,y) 

-2pP,^{x,y)F,b{x,y) 



(53) 



Plsci.x,y) = Fa(s.i){x,y)F(^s,i)c{y,x) 

+ ^Pa{5A) {X, y)P{S,i)c{y, X) , 

PaSci^^y) = Pa{SA}{x,y)Fi^s^,yiy,x) 

- Fa(S,i}{x,y)p{S,z)c{y,x) , 

n^pix.y) = 7^(a,i)(/3j)(a;,|/)7^(/3j)(a,i)(y,x) 

+ ^P(a,j)(0j)(a;,2/)P(/3j)(a,i)(y,a;) , 
F^afji^^y) = P(a,i)(/3,i)(a;,y)^(/3j)(a,»)(y,a;) 

- Fi^a,i){l3,j){x,y)p(i3^j)(a,i){y,x) , 



(54) 



(55) 



(56) 



(57) 



a — (a,«a), b = {l3,ib), and c — (jjic)- On the right- 
hand sides of Eqs. (l5^ -([57 )) . it is summed over indices 
not appearing on the respective left-hand sides. Finally, 
inserting Eqs. ([57)1 and ([SOt into Eq. one finds the 
dynamic equations for F and p: 

rHFB 



dxo-Mfl^ix;G)]F,b{x,y) 



SaJx,z;G)Fcb{z,y) 



yo 



^ac{x,Z]G)pcb{z,y) , 



to 



[iTacdx,-MfJ''{x;G)\pcb{x,y) 



Sac{x,z;G)pcb{z,y) , 



(58) 



(59) 



Vo 



t' t' 

where da; = dxo In deriving these equations, 
the decomposition of G and S into statistical and spec- 
tral parts allows to rewrite time integrals over the CTP 
into simple time integrals. Hence, the right-hand sides 
introduce scattering effects in form of memory integrals 
that render the equations non-Markovian. The form of 
Eqs. (l551) and ([5^ is independent of the order of ap- 
proximation chosen for 72 ■ In the second-order coupling 
approximation introduced above, they form a closed set 
of integro-differential dynamic equations. 

From the dynamic equations (1551) and ([5^ . stan- 
dard quantum kinetic (Boltzmann) equations can be de- 
rived for the mode occupation numbers rta(p,i)- This 
is achieved by a transformation to Wigner space, ne- 
glecting initial-time and non-Markovian effects in a gra- 
dient expansion with respect to the absolute time T = 
{xQ+yo) /2, and making a quasiparticle ansatz cf., e. g., 
Refs. pn[55] for details. The relevance of non-Markovian 
and initial-time effects provided by the full dynamic 
equations has been discussed, for bosonic theories, in 
Refs. ^64.1 57] . For concise discussions of kinetic equations, 
we refer to Refs. 

Following the procedure outlined above, higher-order 
coupling approximations can be derived straightfor- 
wardly. The resulting equations of motion for the two- 
point correlation functions have the form ([55)) and ([5^ 



with modified proper self-energy functions S 
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[G] 




Fig. 4 (Color online) Diagrammatic representation of the 
leading order (LO) and next-to-leading order (NLO) contri- 
butions to /2 in the l/A/" expansion. All statistical factors 
and relative minus signs are omitted. 



3.2 l/Af expansion of the 2PI effective action 



The NLO contribution reads 



rf'^^[G] = -T,[\n[B{x,y-G)]\ (62) 



with 



Bapix, y; G) = SapScix - y) - ^^i7^/3(x, y) (63) 



such that 



Tr[ln[i3(x,2/;G)]] 

iAQ/3 , , 



(64) 



The Lagrangian (U) is symmetric under global U{1) 
transformations of the complex- valued fields ipaM. In 
the special case that all couplings Xap are equaO, the 
Lagrangian has an additional 0{J\f), one has an addi- 
tional 0{N) symmetry in the space of all hyperfine lev- 
els described by the multicomponent field ^a{x). This 
symmetry can be used to derive an expansion of the 2PI 
part /2 of the effective action in powers of the inverse 
number of hyperfine levels N . 

Consider the powers of N in the diagrams in Fig. |4l 
Each vertex contributes a factor of 1 /A/", see Eq. (Hj) . Us- 
ing the 0{Af) symmetry, one can diagonalise the Green 
function for appropriate initial conditions, that is, spin 
balanced mode populations, such that Gap = for 
a ^ p. As a consequence, each loop contributes a factor 
of A/", and all diagrams with one new vertex appearing 
for each new loop contribute to /2 at the same order. 

In the following, we derive the leading-order (LO) 
and next-to-leading-order (NLO) terms of the 2PI 1/M 
expansion. 



r2[G] =r^''[G] + r, 



^NLO 



'G] 



(60) 



The LO contribution is equivalent to one part of the 
Hartree-Fock-Bogoliubov (HFB) diagram. 



rJ^^\G] 



Gaa{x, x)Gbb{x, x) 



(61) 



As there are two sums over a, f3 e {1, ■ . ■ , A/"}, this con- 
tribution is of the same order in Af as the one-loop part 
of the action. Hence, in the limit A/" — >■ oo, the dynamic 
equations contain less terms than in the HFB approxi- 
mation, and the dynamics is entirely mean-field. 



® For the interactions considered here, Eqs. ^ and p2p . 
the fermionic property of the interaction is implemented in 
the fields since 'l^a{x)'I'a{x) = ^a{x)4'a(x) = 0; therefore, 
it is unnecessary to additionally require Xaa = 0, and all 
couplings \ap can indeed be chosen to be equal. 



with 77 defined in Eq. (HH]) . In analogy to Eq. the 
proper self-energy has LO and NLO contributions, 

Sab{x,y,G)^S]^^{x,y-G)+E^t''{x,y-G) + --- , (65) 

with 

E^;^{x,y;G) ^ Sab'-^G,,ix,x)Sc{x ~ y) , (66) 

S^,^''{x,y;G) = -^Aap{x,y-G)Gab{x,y). (67) 

The NLO contribution can be understood as a scattering 
diagram with a resummed vertex 



A 



apix, y) = [Sa-yScix - y) + lla-yix, y, G)) ^ , 



which is defined through the integral equation 

lapix, y;G) = ^ Aajix, z; G)n^fs{z, y) . (69) 

One recovers the Hartree-Fock-Bogoliubov approxima- 
tion by setting the resummed local interaction function 
A equal to the bare coupling, Aa()(x,y) = XapS{x — y)/2. 

The self-energy up to NLO has in general both local 
and nonlocal contributions, the local one being equiva- 
lent to the HFB term, 

Sfrix,y) = Si^ix,y) + S^t''ix,y)\^^^ ■ (70) 

The nonlocal beyond-mean-field contribution to the self- 
energy is given by 



^ab{x,y;G) = Iajix,y;G)'-^Gab{x,y) 



Af 



(71) 



The real functions Mab{x;G), IJ^i^{x,y;G), and 
S'^i,{x,y;G) are all regular in xq and are obtained in 
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terms of statistical and spectral functions as follows: 



and 



Fab {x,x) , 



-:Ia-y{x,y)pab{x,y) 



lbix>y) = -jf-(ja-y{x,y)Pab{x,y) 

+ I^'y{x,y)Fabix,y)j , 



(72) 



(73) 



(74) 



where 



^(/"dz/^,(x,z)i7,^,(z, 



y) 



yo 



dzI^g{x,z)nP{z,y) ] , 



+ W/ dz^:,(x,z)7T^^^(z,y). 



(75) 
(76) 



with and UP defined in Eqs. ([5^ and ([571) . respec- 
tively. As shown in App. [Cl the dynamic equations de- 
rived from the 2PI effective action fulfill crucial conser- 
vation laws including those for the total particle number 
and the total energy. 

In the next section, we focus on a homogeneous gas 
in a box with periodic boundary conditions, in which 
case the equations of motion in momentum space and in 
NLO 1 /Af approximation are given by 

(iTacdxo - Macixo;p)jFcbixo,yo;p) 
dzQ S''^Jxo,zo; p)Fcb{zo, yo; p) 



dzo ^aci^a^ ^o; p)Pcbizo, yo;p), (77) 
(iTacd^o - Macixo;p)jpcb{xQ,ya;p) 
= / dzo S''^^{xo,zo;p)pcb{zo,yo;p) (78) 

J Vn 



with 



Mab{xo;p) = Sab 



p2 A„ 



2m 27V 7k 

^ J Fabixo,xo;k) 



Fddixo,xo;'k) 



(79) 



jjyI^',{xQ,yo;p- k)F„fc(a;o,yo;k) 



I^-y{xo,yQ\p - 'k)pab{xo,yo\'k)) , (80) 



^a6(a;o,yo;p) 



■ J^(l^^{xo, yo;p- k)p„fc(xo, 2/0 ; k) 

+ ^07(2^0, 2/o;p-k)F„f,(xo,2/o;k)^ ■ (81) 
Here, the functions and are given by 

iap{xQ,yo;p) 
= ^n^^p{xo,yn;p) 

'^^oiasixo,zo;p)n^f!{zo,yo;p) 

rVo \ 

-J dzo/f5(xo,zo;p)i7^''^(zo,yo;p)j , (82) 
iaf3i^o,yo;p) 

= ^ 

with 



— ■^^/3(2;o,2/o;p) 



/ dzo/^^(xo,zo;p)i7^''^(zo,2/o;p) , (83) 

Vn 



na0{xo,yo;p) 



(^(a,i)(/3,i)(a;o,yo;P - k)F(^_j)(„_i)(yo,a;o;k) 



rP(a,i)(/3j)(a;o,yo;P - k)p(;3j)(Q,i)(yo,a;o;k)j , 



(84) 



^0/3(2^0, yo;p) 



(p(Q,i)(/3 j)(2^o, yo; P - k)F(0j)(„^i) (yo, a;o; k) 



F(a,i){l3,j) {xo,yo;P- ^)P{l3,j)ia,i) (^0, Xq', k) 



(85) 



where = (27r)"'^ / d'^fc. In one spatial dimension, 
these equations are conveniently written in dimension- 
less variables. This is achieved by defining t = n'lj-it/'m, 
p = p/niD, 7q/3 = rngap/niu, M = mM/nl-^, S = 
m^Z'/n^j) and / = m//n^j-,, where the tilde denotes 
the rescaled quantities, 7q/3 is the dimensionless coupling 
constant and niD denotes the line density. 



4 Nonequilibrium time evolution of a 
one-dimensional Fermi gas 

In this section, we apply the dynamic equations derived 
from the 2PI effective action, in next-to-leading order 
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(NLO) of the 1/M expansion, to the case of a homo- 
geneous ultracold Fermi gas with twofold hyperfine de- 
generacy (denoted as f and \.) in one spatial dimension. 
We study the system in a finite-size box with periodic 
boundary conditions. The Fermi gas is initially assumed 
to be noninteracting and prepared far from equilibrium, 
characterised by a nonequilibrium single-particle mo- 
mentum distribution. The equations of motion, Eqs. (1771) 
and (j78p . are solved in momentum space. We assume the 
interactions to be switched on at the initial time and in- 
vestigate the long-time evolution of the interacting gas 
towards equilibrium. 

The homogeneous one-dimensional gas is taken to 
have a line density tt-id , and the constituents have mass 
m. To identify the relevant combination of parameters, 
it is convenient to rewrite the equations of motion in di- 
mensionless variables, and to introduce the dimension- 
less coupling constant = w^a^/niD. 

Initial conditions 

With initial values for the spectral function p{tQ,tQ;p) 
and the statistical propagator F{to,to;p), Eqs. (iTTl) and 
([75]) describe the time evolution of the two-time correla- 
tion functions including the momentum distribution 

n„(<,p) = i(l-f^(„^,)(„^,)(i,t;p)). (86) 

In the following, we will choose different initial momen- 
tum distributions na{tQ^p) away from thermal equilib- 
rium. Furthermore, we choose the initial coherence be- 
tween different spins as well as the initial pair correlation 
function to vanish, 

(#t(x,to)<^^(x,to)) =0 fora^/3, (87) 

(#Jx,to)'^/3(x,to)) =0. (88) 

For a = /3, the equal-time pair correlation function al- 
ways vanishes, which is in accordance with the conser- 
vation of the total particle number and a direct conse- 
quence of the equal-time property of the spectral func- 
tion (see App. [B|) . For a 7^ /3, a nonzero initial pair cor- 
relation function would account for BCS-type pairs and 
imply a nonzero variance of the total particle number. 
The above initial conditions require 

F{a,ia){i3,it){to,ta;p) ^0 ioriaj^ib- (89) 

Combining Eqs. ((55)) and yields the initial con- 
dition 

= ^-^a(io,p). (90) 

The equal-time property of the spectral function, 
Eq. (|112l) . requires 

Pab{to,to;p) ^^Tab- (91) 



For a homogeneous gas and a (p — p)-symmetric 
initial state, F and p are invariant under p — !■ — p at all 
later times. 

4.2 Numerical method 

We numerically solve the equations of motion (j77p and 
(|78p together with the nonperturbative integral equa- 
tions for the self-energies (|80)) and (|8T|) . Due to the mem- 
ory integrals, computations are costly, and in computing 
the results shown below, we kept a finite memory kernel 
at longer times, checking that an increase in the mem- 
ory time did not change the results. We also carefully 
checked that a change in the size of the box and the 
number of the momentum modes does not lead to a sig- 
nificant change in the results presented in the follow- 
ing, except for finite-size (infrared-cutoff) effects taken 
into account explicitly. For example, we chose the in- 
teraction strength sufficiently weak such that the oc- 
cupation numbers of modes close to the cutoff defined 
by the grid size are small enough not to give rise to 
ultraviolet-cutoff-dependent effects. For the results pre- 
sented here, we discretised the momentum space on a 
lattice with Ng = 128 sites and periodic boundary con- 
ditions. Thus, the momentum modes corresponding to 
the lattice Laplacian are Pj/niD = 2NsSm{jTr/Ng)/N, 
where j € {-Ns/2 + 1, ~N,/2 + 2, . . . , Ns/2} and N = 
j "■a(Pj") is the total particle number. 
The diagonal time steps were implemented according 

to 

F{tn+i,tn+i;p) ^ F{tn,tn;p) 

+ (^F{tn+l,tn;p) - F{tn,tn;p)j 

+ (F{tn,tn+i;p) ~ F{tn,tn]p)) 

(92) 

in order to ensure that the equations of motion, which 
are implemented on discrete grid also in the time do- 
main, satisfy the same global U{1) symmetry as the con- 
tinuous equations. In this way, the total particle number 
is numerically exactly conserved. 

4-3 Results 

To begin with, we studied the time-evolution for two ini- 
tial states with the same total particle number and total 
energy (and interaction strength 7 = 4), but different 
far-from-equilibrium momentum distributions. The two 
initial momentum distributions are shown in Fig. [5^. 

Once the system evolves in time, multiple scatter- 
ing events lead to a redistribution of momenta until the 
system reaches an equilibrated state. For runs A and B, 
Fig. [5}) shows the momentum- mode occupation numbers 
as a function of time t for six of the momentum modes. 
We find that the equilibrated states at late times have 
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the same momentum-mode occupation numbers. Note 
that a conservative estimate of the recurrence time on 
the basis of the slowest osciUating discrete mode of the 
free gas is two orders of magnitude larger than the shown 
total evolution time. The final {tn'f^/m ~ 10) momen- 
tum distributions are also shown in Fig. [5^. However, 
to see that the final states are absolutely identical, one 
has not only to look at the momentum-mode occupa- 
tion numbers, but also to compare the nonlocal-in-time 
behaviour of the two-point functions Fabixo,yo]p) f^nd 
Pab{xo,yo]p) of the two runs for same center time co- 
ordinates X = (xq -|- yo)l2 but different relative time 
coordinates s = — j/o- Exemplarily, i)(-f- s;p) 
at late times {n\j^X/m = 10) is depicted in Fig. [S}: for 
four of the momentum modes of runs A and B. The data 
for the two runs lie on top of each other; thus, the final 
correlation functions of runs A and B are indeed identi- 
cal. 

Furthermore, to observe that the final momentum 
distribution has the form of a Fermi-Dirac distribution 
within the range of momenta considered, it is convenient 
to look at the inverse slope function a = ln[l/n^ — 1] 
rather than the occupation numbers directly. When sub- 
stituting the inverse slope function for the exponent of 
the Fermi-Dirac distribution, i.e., 



(a) 



1 



exp{a-[tj(p)]} -f- 1 



(93) 



it is apparent that cr, as a function of mode energy u;(j»), 
reduces to a straight line when the occupation num- 
ber n-^{t,p) is a Fermi-Dirac distribution. Fig. [5^ shows 
the inverse slope function as a function of mode energy 
uj{p) = for run B at various times t. At late times, 
where the inverse slope function becomes a straight line, 
one can extract the inverse temperature j3 and the chem- 
ical potential /i according to cr(uj) = /3[uj — fi). 

A thermal state, however, does not necessarily re- 
quire that the single-particle momentum distribution co- 
incides with a Fermi-Dirac distribution. For the state 
to be described as a (grand-)canonical ensemble, a 
necessary condition is the Callen-Welton fiuctuation- 
dissipation relation |67il68l 



FaaiXo;uj,p) 



-i[l/2~nFu{^-'H)]paa{XQ;LL>,p) (94) 



which connects the statistical propagator F and 
the spectral function p. Here, Faa{X;Lu,p) ~ 
JdseMiLjs)Fo,^{X + s/2,X - s/2;p), F^^{t,t'-p) = 
([#^(t,p),#a(t',p)]_)/2, and similar for paa{t,t';p) = 
i{[Wl{t,p),iFait',p)]+). We emphasise that Eq. ^ is 
valid for a grand canonical state irrespective of whether 
the system is interacting, i.e., whether single-particle 
modes are eigenmodes of the Hamiltonian, or not. 

Figure [Sb shows the emergence of the fiuctuation- 
dissipation relation at late times for run B, where we 
introduced the fraction 
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Fig. 5 (color onhne) Comparison between run A and run 
B. (a) Initial (nf^f/m — 0) and final (n\-Qt/m = 10) 
momentum-mode occupation numbers n-|-(t,p) = n^(t, —p) = 
ni{t,p). The initial states are set up such that both runs have 
the same total particle number and total energy. Therefore, 
the equilibrated states at late times have the same momen- 
tum distribution, (b) Numerically determined momentum- 
mode occupation numbers as a function of time t for some 
of the momentum modes. Note that time and occupation 
numbers are shown on a logarithmic scale. Note that a 
conservative estimate of the recurrence time on the basis 
of the slowest oscillating discrete mode of the free gas is 
two orders of magnitude larger than the shown total evo- 
lution time, (c) Nonlocal-in-time behaviour of the enve- 
lope functions s;p) of the statistical propaga- 
tor %,i)(^,i)(X, s;p) = F(^,i)(^,i)(X, s;p)cos[p^s/(\/2m)] at 
late center times (nf^X/m = 10) for some of the momentum 
modes. The data lie on top of each other, which proofs that 
the final states of both runs are identical. 



/ = \F^^{Xq]uj,p)/ pA;-,;{Xa;uj,p) 



as a function of the frequency w. / is shown in a region 
around the peaks of the statistical and spectral functions 
where the argument of the logarithm ln[(l/2 — /)~^ — 1] 
is positive. Outside this region it oscillates around zero 
due to finite evolution time after the quench. Propagat- 
(95) ing the dynamic equations further reduces these oscil- 



Far-from-equilibrium dynamics of an ultracold Fermi gas 



13 



Si, 



(b) 



I 

o 



' 1 


' 1 ' 1 


'I'll' ^ 






<^^\ = o • 






i = 16 ' - 






i = 32 ' 






i = 48 • - 


, 1 


1 , T 1 


i = 64 ■ 

I.I. 



1^ 




Fig. 6 (color online) Further analysis of run B. (a) Fermion 
distribution nf{t,p) as a function of mode energy uj{p) in 
units of the Fermi energy ljf = uj{pf) at various times t — 
n1j^t/m. Plotted is the inverse slope function a — ln[l/n-f — 1], 
which reduces to a straight line when nf{t,p) approaches a 
Fermi-Dirac distribution, (b) Fluctuation-dissipation relation 
at late times (nlj^t/m — 10) for several momentum modes 
Pi. The fraction / is defined in Eq. H95p . (c) Time evolution 
of the potential and total energy. The reduction of the initial 
potential energy is due to the built-up of additional kinetic 
energy. 



lations. Hence, according to the fluctuation-dissipation 
theorem, the system is approximately thermahsed over 
the depicted range of energies. 

The redistribution of the initial kinetic and potential 
energies during the time evolution is shown in Fig. |6]:. 
Since the total energy is calculated by summing the nu- 
merically determined kinetic and potential energies, this 
plot also highlights the numerical accuracy of the con- 



servation of the total energy. We have explicitly checked 
that extending the size of the memory kernel further 
does not change our results. 

In conlusion, the chosen initial conditions for runs A 
and B allovi' the single-particle momentum distributions 
to thermally equilibrate over the considered range of mo- 
menta. We find that the final state is determined by the 
values of the conserved quantities in the initial state. All 
other information about the details of the initial state is 
lost during the evolution. 

Next, we investigate a few thermodynamic properties 
of the equilibrated interacting Fermi gas. For this pur- 
pose, we performed additional runs with different initial 
energies but the same total particle number and inter- 
action strength as in runs A and B. Exemplarily, the 
initial momentum distributions of two of these runs are 
shown in Fig. [7^ together with the one from run B. Even 
though all of the runs virtually reach a stationary state 
(Fig. [7}d) within the times that are numerically accessi- 
ble before the time discretisation leads to a break-down 
of the energy conservation, the runs with lower initial 
energies do not fully settle to a Fermi-Dirac distribu- 
tion. As depicted in Fig.[7j;, the lower momentum modes 
reach a Fermi-Dirac distribution, but the higher momen- 
tum modes still show an excess population. Keeping the 
total particle number constant, a further reduction of 
the population in the higher momentum modes will not 
significantly alter the population in the lower momen- 
tum modes since the former are already populated much 
less than the latter. Therefore, we can extract tempera- 
tures and chemical potentials from a fit of the lowest 14 
momentum modes to a Fermi-Dirac distribution as it is 
shown in Fig. [T]:. 

The so found temperature dependence of the late- 
time kinetic energy i?^in^ heat capacity C'v = 
ksdEtot/ dl5~^ at constant volume, and the chemical po- 
tential /J, are shown in Fig. [5] At high temperatures, the 
quantities of the interacting gas converge towards the 
results for an ideal Fermi gas. However, at low temper- 
atures, they significantly deviate from those of an ideal 
gas due to the finite coupling constant as shown in the 
insets. Note also that the results are sensitive to the fi- 
nite size of our system. Therefore, also the results for the 
ideal Fermi gas in discrete momentum space differ from 
those of an ideal Fermi gas in the thermodynamic limit, 
where the momentum space is continuous. 

Finally, we investigate the overpopulation of the oc- 
cupation numbers for runs C and D at high momenta 
as compared to a Fermi-Dirac distribution. As depicted 
in the inset of Fig. [Tj:, the tail in the momentum dis- 
tribution is characterised by a power-law n^(p) on p~'^ 
with K ~ 4.4. This exponent does not change when dia- 
grams of order A'^ are included in the effective action, see 
Fig. |9l Note, however, that it is crucial to go beyond the 
LO l/JV and the HFB approximations of the effective 
action since both of them exclude multiple scattering 
events and thereby leave all momentum-mode occupa- 
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Fig. 7 (color onhne) (a) Initial momentum distributions 
of runs B, C, and D. The runs have the same total particle 
number but different total energies, (b) Time evolution of 
the occupation numbers for some of the momentum modes 
of run D. Similar to all other performed runs, a stationary 
state is reached at late times, (c) Inverse slope functions for 
runs B, C, and D at late times. The black dashed lines are 
fits to the lowest 14 momentum modes from which the tem- 
peratures and chemical potentials are extracted. The inset 
depicts the power-law tail in the momentum distributions 
for low energies. 



tion numbers unchagend during the time evolution, of. 
Eq. (HH) and the ensuing discussion. 

As pointed out before, this power-law tail suggests 
but does not prove that the Fermi gas approaches a 
non-thermal state. Therefore, we have a closer look 
at the fluctuation-dissipation relation, Eq. The 
lower panel of Fig. [TU] shows the inverse-slope function 
ln[(l/2-/)-i-l] of the fraction / defined in Eq. ^ for 
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Fig. 8 (color online) Temperature dependence of (a) the 
mean kinetic energy per particle, (b) the heat capacity per 
particle, and (c) the chemical potential. The results for the 
interacting Fermi gas from the 2PI 1/M NLO runs are shown 
as red dots. For comparison, the exact results for the ideal 
Fermi gas are shown both for a continous and a discrete mo- 
mentum space. The differences between the interacting Fermi 
gas and the ideal Fermi gas in discrete momentum space are 
shown in the insets. 



run C, the lower panel of Fig. llll for run D. / is shown in 
a region where the argument of the logarithm is positive 
~ outside this region, it oscillates around zero as a result 
of the finite total evolution time after the quench. 

In run C, the inverse-slope function is a straight line 
over the region of relevant uj and therefore corresponds to 
a Fermi-Dirac function. As in run B, the system is ther- 
malised over the depicted range of energies, in spite of 
the signs of a power- law tail in run C. This can be under- 
stood by considering the spectral function in the upper 
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Fig. 9 (color online) Momentum distributions at niY)t/m = 
6.4 for runs with the same initial state as run D but in dif- 
ferent approximations of the effective action. 



Fig. 11 (color online) Same as Fig. [10] but for run D. In 

run D, the system does not thermalise. Note the exponential 
decay of the spectral functions away from the peaks. 
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a Fermi gas in the BCS theory, a p power-law tail of 
the single-particle momentum distribution. 

However, run D performed at even lower energy 
shows that the system does in general not thermalise to 
a grand-canonical ensemble: even though the momentum 
overpopulation is again largely produced by the contri- 
butions from the Fermi sea, see Fig. [TT] (upper panel) , 
also the fraction / shows a power-law tail ~ vio- 
lating the fluctuation-dissipation theorem, see Fig. [TT] 
(lower panel). Despite this, the equilibrated momentum 
distribution oc p~^'^ at low total energies is still mainly 
due to the second peak in the spectral function while the 
contribution from the nonthermal power-law tail of / is 
suppressed by another 4.5 powers of p. 



Fig. 10 (color online) Top: Spectral functions as a function 
of frequency at late time Xo = 18.9 n^^m in runs C for five 
of the momentum modes p^. Bottom: Inverse-slope function 
of fractions / of the statistical F divided by the spectral 
function p at Xo — 18.9 n^^m, for the same five momentum 
modes. Black lines indicate Fermi-Dirac distributions with /3 
and jj, as in Fig. [7];. 

panel of Fig. [TU] the second peak at negative frequencies 
picks up extra contributions from the Fermi sea thereby 
causing the power-law overpopulation at high momenta. 
Although the area under the negative-cj peak is reduced 
by a factor of ^ 10^'*, it is multiplied by 1 on the filled- 
sea side of the Fermi-Dirac function while the positive-w 
peak multiplies the exponentially suppressed tail of the 
Fermi-Dirac function. Thus, one may preconclude from 
run C that the system thermalises to a grand-canonical 
ensemble, with the eigenmodes of the strongly interact- 
ing system at low temperatures being superpositions of 
particles and holes. This contains reminiscence of the 
Bogoliubov depletion at zero temperature that gives, for 



5 Conclusions 

We have summarised the description of far-from- 
equilibrium dynamics of ultracold Fermi gases in 
terms of Kadanoff-Baym equations for two-point many- 
body Green functions derived from the two-particle- 
irreducible (2PI) effective action in nonperturbative ap- 
proximation. The approach allows to handle both mean- 
field and beyond mean-field approximations on the same 
footings and to derive approximations far beyond mean- 
field and perturbative kinetic approaches in an elegant 
way. Obtaining approximations on the level of the 2PI 
effective action ensures the conservation of energy irre- 
spective of the chosen truncation as well as other vital 
quantities as the total particle number. This forms a pre- 
condition for the applicability of the approach for long- 
term evolution and equilibration. Implicitly contained 
higher-order correlations render the dynamic equations 
nonlocal in time, causing a non-Markovian scattering in- 
tegral. 
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Beyond the mean-field truncation wliicli excludes ef- 
fects of scattering between quasiparticles, we consider 
the nonperturbative expansion in inverse powers of the 
number of internal or spin degrees of freedom A/". This 
is possible when the external potential and collisional 
interactions show an A/"- fold spin degeneracy, i. e., if the 
Hamiltonian is symmetric under the orthogonal group 
0{Af). We have considered a system of Fermions inter- 
acting through s-wave collisions between different spin 
components. 

As a specific example, we have studied the long-time 
evolution of a homogeneous, one-dimensional, twofold 
spin-degenerate Fermi gas with an initial far-from- 
equilibrium momentum distribution. We have numeri- 
cally solved the dynamical equations in next-to-leading 
order in the 1 /Af expansion of the 2PI effective action. 
Results presented for this extend upon the work pre- 
sented in [43) . They give that, within the truncation 
considered, the one-dimensional gas dephases and equi- 
librates following an interaction quench. For sufficiently 
high total energies, the system is found to equilibrate 
to a state with thermodynamic properties like chemi- 
cal potential and specific heat given by those of a ther- 
mal ideal Fermi gas. In contrast, close to the Fermi en- 
ergy the equilibration leads to non-thermal power-law 
momentum distributions pointing to the appearance of 
many-body quasiparticle modes. Furthermore, we found 
a violation of the fluctuation-dissipation relation for a 
grand-canonical ensemble which gives a strong indica- 
tion of a nonthermal equilibrium state. A future task is 
to extend our study to account for the limit of low en- 
ergies near zero temperature, in particular for features 
of the Tomonaga-Luttinger low-energy fixed point, in or- 
der to clarify the transition regime between high and low 
energies. 

Acknowledgements The authors would like to thank 
C. Bodet, M. Holland, S. Jochim, S. Kehrein, V. Meden, 
B. Nowak, J. M. Pawlowski, A.M. Rey, and D. Sexty for in- 
spiring and useful discussions, and JILA and the University 
of Colorado for their hospitality. They acknowledge the sup- 
port by the Deutsche Forschungsgemeinschaft, as well as the 
support of the Alliance Program of the Helmholtz Associa- 
tion (HA216/EMMI). M. K. thanks the Heidelberg Graduate 
School for Fundamental Physics (HGSFP) and the German 
Academic Exchange Service (DAAD) for financial support 
during his stay at JILA. 



A Grafimann variables 

In this appendix, we outline properties of Grafimann 
variables that are needed in the context of our discus- 
sion of nonrelativistic fermionic path integrals and not 
commonly found in textbooks. For a more detailed dis- 
cussion of Grafimann variables in the context of path 
integrals, we refer, e.g., to Ref. [69j. 



A set of complex Grafimann variables 9i, i e 
{1, . . . , n}, satisfies 

e,e, + eje, = o, i,j e {i,2,...,n}. (96) 

This implies 

0^9, = , (97) 

with no summation over i. 

For any pair of Grafimann variables 6 and z?, complex 
conjugation (denoted by an asterisk) is defined by 

(dd)* = ■&*9* , (98) 

which ensures the reality condition {09*)* — 90*. 

The real and imaginary parts of a complex Grafimann 
variable 9 are defined as 

i=0i=Re[0]^l(^? + r), (99) 

±0,=lni[0]^^^[0-9*). (100) 

From this definition of 9i and 02, and property ((96l) of 9 
and 9* follows immediately 

9101= 0202 = 0, (101) 

10102 = -i029i= 9*9, (102) 

which implies that 1611 6*2 rather than 0i02 is real. This 
motivates to introduce the bared notation, 

01 = -102 , 92 = i0i , (103) 

so that 

9,0 J = -0J0, e M , i,je {1,2}, (104) 

and 

0*0 ='^{0101 + 0292)- (105) 
B Two-point Green functions 

In this appendix, we would like to summarise symmetry 
properties of the two-point Green function G used in our 
derivations, as well as the statistical and spectral com- 
ponents of the two point function, F and p, respectively. 

The two-point Green function G is defined as the 
time-ordered expectation value of two fields at two 
points in space-time, 

Gab{x,y) = {TcMxWbiy)) , (106) 

where 7c denotes time-ordering along the closed time 
path C. Since there are no fermionic field expectation val- 
ues, G is automatically connected. For Grafimann fields, 
the time-ordering is defined as 




M=^)^b{y) if sgnc(a:o - J/o) - 1 (^07) 
-Wi,{y)'ifa{x) if sgnc(a;o - Vq) = -1 , 
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where sgn^^xo ~ yo) denotes the sign function along the 
time path C and evahiates to 1 (—1) if xq is posterior 
(prior) to yo. 

The path integral formulation naturally accounts for 
time-ordering in expectation values of operator prod- 
ucts. At equal times, however, the time-ordering is ill- 
defined. With the above definition the two-point Green 
function is singular at a;o = yo - For numerical implemen- 
tations, it is convenient to make the singularity explicit 
and to decompose G into its non-singular spectral and 
statistical components. The spectral function p contains 
information about the spectrum of the theory, i.e., the 
energies and decay times of the states, and the statisti- 
cal propagator F accounts for the respective occupation 
numbers. They are defined as 



Fabix^y) 
Pab{x,y) 



i([#,(x),#,(y)] , 



where 



denotes the commutator and 



(108) 
(109) 
the an- 



ticommutator. The decomposition identity then reads 



C Conservation laws 

For a theory of dynamics to be physically meaningful, 
it is crucial to respect the conservation laws prescribed 
by the symmetries present in nature. For example, for 
a system that is closed with respect to the exchange of 
particles and energy with its surroundings, a theoretical 
description of its dynamics needs to refiect the conser- 
vation of particle number and energy irrespective of the 
approximation chosen. Moreover, if the system relaxes 
to an equilibrated state in the long-time limit then this 
state is determined by the values of the conserved quan- 
tities in the initial state. In the following, we show that 
the total particle number and energy are conserved at 
any order of truncation of the two-particle irrecducible 
(2PI) effective action. To take into account other inde- 
pendent conserved quantities, a generalization of the 2PI 
effective action approach to nPI effective actions is re- 
quired [72] . 



C.l Particle number conservation 



Gab{x,y) ^ Fab{x,y)~-pab{x,y)sg-nc{xo - yo) . (110) 

The decomposition of G into F and p has the advantage 
that the time ordering, and thus the discontinuity at 
Xo = yoi is accounted for explicitly by the sign function. 

We close this appendix with some useful symmetry 
properties of the two-point functions: 

Gab{x,y) = TbcGcd{y,x)Tda = {-'^)''"'^'"'Gi^{y,x) , 

(111a) 

Fab{x,y) =ncFcd{y,x)Tda =(-l)'°+"'Fba(j/,a;), 

(111b) 

Pab{x,y) = ncPcd{y,x)Tda = {-iy''^'-''Pla{y,x) , 

(lllc) 

where a — {a, 3 ~ ia) and ia — 3 — ia- Note that in the 
definitions of p and F the operator product is not time 
ordered and cannot be directly evaluated in the path in- 
tegral formalism. However, expectation values of equal- 
time anticommutators of Grafimann variables can be 
evaluated using the Bjorken- Johnson-Low theorem [701 
[71] . For the spectral density function, one findfQ: 

Pab{x,y)l^^^^^\T,b5^^\^-y). (112) 

In the derivation of the equations of motion we made use 
of this anticommutation relation in the identity 



Pab{x, y)d.^a sgnc(a:o - yo) 



'2'Pab{x,y)5c{xo 
2iTab5c{x - y) ■ 



yo) 



(113) 



The anticommutation relations are consistent with the 
anticommutation relations in the canonical quantization ap- 
proach. 



According to Noether's theorem, the symmetry of the 
Lagrangian (j?]) with respect to global U{1) transforma- 
tions of the complex fields ^pa implies the conservation of 
the total particle number. This can be shown as follows. 

From the stationary condition (j27L one can con- 
struct the vanishing expression 

= -2ir,, / -ffiLG,a(y, x) . (114) 



"'Jy 5G,b{y,x) 

Substituting the right-hand side of Eq. (|^ for F into 
this expression, one finds 

5.o^(^)+V.j(x) = -2ir,, / ,^^"f^\ GUy.x) (115) 



with 



SGcbiy,x) 



nnt[G]=F[G] + ^Tr[Go'G] 



n{x) = y(i ^ Gaa{x,x)j 



i{x) = i7^V;,GtJz,y) 
2m 



(116) 
(117) 
(118) 



z=y=x 



where sums over a = {a,ia), b = (/3, «&), and c = 
(7,ic) are implied. Like the underlying classical action 
S, Eq. ([U, the effective action, and hence also /Int, is 
a singlet under 0{M) rotations and parameterized by 
the field G. All functions that are singlet under 0{N) 
rotations can be built from the irreducible, i.e., in spin- 
index not factorizable, invariants tr[G"] with n < Af. 
Here, tr[-] applies to the spin indices and does not in- 
clude an integration over space-time coordinates, e.g., 
tr[G^] = Gab{x,y)Gbc{y,z)Gcaiz,x). Thus, the inte- 
grand of Eq. (jllSp is symmetric under an exchange of a 
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and p. Tab is antisymmetric under an exchange of a and 
p. Hence, summing over a and b sets the right-hand side 
of Eq. (jll5l) to zero, and Eq. (IllSp becomes the continu- 
ity equation. These symmetry considerations are inde- 
pendent of the specific approximation to -Tint- Thus, for 
a closed system, the total particle number is conserved 
at any order of truncation of the 2PI effective action. 



C.2 Energy conservation 

For time independent interactions and external poten- 
tial, the Lagrangian is time translation invariant. This 
implies energy conservation. Here, we consider the in- 
variance under general translations in continuous space 



and time that vanish at the boundary, 



-e''(x). 



where e'' (x) is a time- and space-dependent infinitesimal 
{d + l)-vector. To leading order in e, the Green func- 
tion transforms under these translations as Gab{x, y) — 
Gab{.x,y) + e''ix)d^Gab{x,y) + e'' {y)dyGab{x,y), where 
9^ = d/dx'^. One can show that under these transfor- 
mations the variation of the 2PI effective action F can 
be written as r[G] ^ r[G] + Sr[G], with 



Sr[G]= f T^''{x)dleAx). 

J X 



(119) 



Since, by virtue of the stationarity condition (l?7|) . the 
variation SF vanishes for all solutions of the equation 
of motion for G, an integration by parts shows that 
T^^ is the conserved Noether current for the time-space- 
translations: 



5r[G] 



e,ix)d^^T^"ix)^0. 



(120) 



T'^'^(x) is identified as the energy-momentum tensor, 
and the conservation law for total energy is expressed 
as d^Tf^^ix) =0oi dtj S'^hT°°{t,-si) = 0. 

As sketched in the following, the energy momentum 
tensor can be calculated by standard techniques for any 
truncation of the 2PI effective action. The energy density 
is found to be 



Too(x) = £kin(a;) +epot(a;) 



with 



(121) 



(122) 



-pot 



(^) = ^rijx){l n„{x)) ll^Jx,x) , (123) 



ekin(a;) = / S{x - y)H„Jx)n„Jx,y) , 

K 

2J\f' 

where sums over a and /3 are implied. The conserved 
total energy is then given as the spatial integral over the 
above density. 

We close this appendix with a sketch of the deriva- 
tion of the energy momentum tensor for an ultracold 
Fermi gas in d spatial dimensions, described by the 2PI 
effective action in next-to-leading order (NLO) 1/Af ap- 
proximation. To this end, the variation of the effective 



action F under space-time translations is split into one- 
loop and higher order terms, 

5F[G] ^-^S{Tr[\nG-^ +G^^G])+6F2[G]. (124) 

To obtain the contribution to T^"^ arising from the one- 
loop term, we use Eq. pi9D and observe that the term 
^Tr [inG-i] in Eq. (|TM1) does not contribute to T'^'". 
The term Tr [Gq ^G] can be written as 



s(-^Tr[G^\x,y)G{y,x)] 



-7; / Go.^b(a:,y)5Gf,„(y,x) 

^ Jxy 



(125) 



with (the potential is assumed to vanish, Vo^t.ap = 0) 
iG'o,L(^. y) = 5{x y) {iTabdl + Sab^dldl) , (126) 



SG 



a,{x,y) = [e\x)d^,+s\y)dl)Ga,{x,y). (127) 



Integration by parts and the identity J^dl[6{x — 
y)G{y,x)] = J^S{x — y){d^+dl)G{y,x) allows to rewrite 
the r. h. s. of E^q. UTI^i in the form of Eq. ([TTOl) . The re- 
sulting contribution to the energy density Too{x) is 

e^inix) S{x - y)5abH]^p{x)G,a{x, y) . (128) 

To obtain the contribution from the term F2 to the 
energy-momentum tensor, it is more convenient to use 
a method known from field theory on curved space- 
time, where, for a space-time-dependent metric g^'^ , the 
energy-momentum tensor is given by |73| 



Tf_,y{x) 



2 5F[G, 



(129) 



with ^y—g{x) denoting the square root of minus 
the determinant of g^^(a;). We apply this to F2 
in NLO of the 1/Af e xpansion, cf. Eqs. ([SOHS^ . 
with J = J dxy^—g{x). Making use of the identity 
Sy/^^/Sg^'^iy) = ^^^^g^,{x)S{x - y)/2 in flat 
space with g^v{x) = diag{l, —1, —1, —1}, one finds 

~ .9p.(^)^^aa(^, X)G,,{X, X) , 

(130) 

— - — ^ , , = —-qi,Jx)Iaa{x,x) , (131) 

where the function I^p{x,y) is given by Eq. (IMl) . Fi- 
nally, using Gaa{x,x) = Faa(x,x) = 1 — 2na(x) allows 
to rewrite Eqs. (IT^ . ([T^ and to obtain the re- 

sults (|122l) and (|123p . where also the vacuum energies 
are removed. 
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